function [ F ] = evaluateInteriorDeformationGradient( x, y, z, mesh, theta )
%EVALUATEINTERIORDEFORMATIONGRADIENT Summary of this function goes here
%   Detailed explanation goes here

n = size(x, 1);
nTri = size(mesh.f, 1);

F = zeros(3, 3, n);

for i = 1 : n
    
    p = [x(i), y(i), z(i)];
    
    for j = 1 : nTri        
        
        q0 = mesh.v(mesh.f(j, 1), :);
        q1 = mesh.v(mesh.f(j, 2), :);
        q2 = mesh.v(mesh.f(j, 3), :);
        
        T = triarea(q0, q1, q2);
              
%         qn = mesh.fn(j, :);
        
        v0 = getKernelGradient(p, q0);
        v1 = getKernelGradient(p, q1);
        v2 = getKernelGradient(p, q2);
        
        Gu = tensorVectorContraction((v0 + v1 + v2) * T / 3.0, theta(3*j-2:3*j)');
        
%         temp = tensorVectorContraction(v2, theta(3*j-2:3*j)');
        
%         temp(:,2)';
        
        F(:, :, i) = F(:, :, i) + Gu + Gu';
 
    end
    
end

% F = -1 * F;

end

